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Description 


Just  about  the  most  common  multi-server,  batch-arrival  queueing  system  encountered  in  practice  is  the 
bulk  Poisson  input  and  constant  service  model,  symbolized  by  Mx / D/c.  A  prototypical  application  arises 
in  vehicle  loading  problems:  batches  arrive  as  a  Poisson  process  to  be  loaded  onto  numerous  vechicles  at  a 
fixed  rate.  The  analytic  development  of  the  model  for  the  bulk-arrival,  multiserver  queue  Mx / D/c  is  given 
in  Chaudhry  and  Templeton  (1983).  Further  computational  aspects  of  the  computation  of  the  distribution 
of  customers  have  been  discussed  in  Chaudhry  et  al.  (1988).  Smith  (1987)  considered  only  an  approximate 
solution  to  a  special  case  of  Mx  j  D/c,  viz.,  the  non-bulk  M/ D/c  queue.  We  present  here  an  exact  numerical 
solution  for  the  distribution  of  the  number  in  the  system  for  the  queueing  model  Mx / D/c. 

Many  applications  of  the  model  Mx  f  D/c  may  involve  high  values  of  the  model  parameters,  particularly 
if  the  model  is  applied  to  large-scale  systems  in  telecommunication  and  transportation.  The  purpose  of  this 
paper,  therefore,  is  to  complement  the  work  reported  by  Chaudhry  et  al.  (1988)  by  producing  a  refined 
version  of  the  program  used  for  that  paper.  The  program  will  run  for  high  values  of  the  bulk  arrival 
parameters  (group  sizes  <  100),  high  and  low  values  of  c  and  p  (1  <  c  <  100,  .01  <  p  <  .9).  Also,  limited 
testing  has  indicated  that  the  program  should  work  for  more  extreme  values  as  well. 
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The  probability  generating  function  (p.g.f.)  for  the  bulk-arrival,  multi-server  system  Mx / D/c  is  in¬ 
verted  nnmerically  by  finding  the  roots  of  the  characteristic  equation  (c.e.)  of  the  p.g.f..  The  roots  are  found 
here  using  a  special  algorithm  that  combines  the  bisection  and  Muller’s  methods.  The  probabilities  for  the 
number  in  the  system  are  calculated  by  using  recursive  formulas  given  by  equations  (2), (3), (6), (7),  and  (8). 
Then  if  needed  (for  example,  when  c  is  “large*  or  there  is  a  ‘large*  number  of  expected  arrivals  during 
service),  the  formulas  given  by  (4), and  (5)  can  be  used  in  an  iterative  way  to  increase  the  accuracy  of  the 
probabilities. 

Formulas: 

Using  the  notation  found  for  example  in  Chaudhry  and  Templeton  (1983): 

A  =  arrival  rate;  b  =  duration  of  service  time  (constant);  c  =  number  of  servers; 
ay  =  Prob(  group  size,  X  =  j);  3  =  mean  batch  size  =  £  j ay  ;  p  =  utilisation  factor  =  Adfc/c; 
xn(t)  =  Prob[  total  number  of  arrivals  during  (0,  t)  =  n];  Pn  =  Prob[  n  customers  are  in  the  system] 
The  underlying  c.e.  of  the  model  is: 

-  exp{6A(A(z)  -  1)}  =  0  (1) 

which  has  exactly  (c  —  1)  distinct  roots  within  and  one  root  on  the  unit  circle,  |z|  =  1.  Let  the  roots  be 
Zi , ...,  Zg—i,  Ze  —  1.  Chaudhry  et  al.  (1988)  show  that  the  first  c  steady-state  system  size  probabilities  in 
Mx/D/c  are  given  by 


p«=cd-p)n(':  , 

»= 1  v  7 

(2) 

Pn  ~  bnPo  y  1  ^  fl  <  C 

(3) 

with  bn  —  coefficient  of  zn  in  f]*=l(l  “  zt zi)-  Then,  using  the  closed-form  equations  for  icn  (also  given  in 
Chaudhry  and  Templeton,  1983)  and  the  steady-state  Chapman-Kolmogorov  (C-K)  equations  written  as 

fl>=Mi)£Pm,  (4) 

m=0 

c  c+n 

Pn  —  ^nfi)  ^  Pm  +  ^  ^  Pm^n-m+e{b)t  f\  —  1,2, ...  (5) 

m=0  m=e+l 
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we  find  that 


Pe  =  P0exp(c/>/a)  -  ]T]  Pm, 

m=0 

p  _  ~  Tl(fe)  Em=0  -P m 

C+l  *o(*) 


r,  Pn~  Xf\(b)  Z)m=0  ~  £m=c+\  Pm*n-m+c(b) _ „  , 

Pe+„ - - ,  n  —  2,3,. 


Method  for  Calculating  Pn 

The  probabilities  are  calculated  by  first  finding  the  roots  of  (1)  and  then  using  the  recursive  relations 
given  by  equations  (2), (3), (6), (7),  and  (8).  The  calculations  proceed  until  either  the  sum  of  the  probabilities 
in  the  vector  {P„}  is  very  close  to  one  or  the  maximum  dimension  of  the  probability  vector  {Pn}  i»  reached. 

The  C-K  equations  (4)  and  (5)  are  applied  to  the  probability  vector  as  an  iteration  in  order  to  increase 
the  accuracy  of  the  probabilities  whenever  the  recursive  relations  fail  to  achieve  satisfactory  results.  The 
probability  vector  is  then  iterated  until  either  one  of  two  conditions  is  met.  If  we  let  { PQBn }  =  probability 
vector  before  an  iteration,  and  {Pn}  =  probability  vector  after  an  iteration,  then  the  basic  condition  is: 

me 

DIFF  =  £  \Pn  -  PQBn |  <  DIFM 

n=0 

where  size  is  the  sise  of  the  probability  vector  {Pn},  and  DIFM  is  a  tolerance  to  be  given  by  the  user  as 
an  input  parameter.  The  other  condition  is  that  the  number  of  iterations  cannot  exceed  the  input  parameter 
MAXIT. 

In  most  cases,  the  recursive  method  will  produce  very  accurate  results  and  the  iterative  method  will 
not  be  needed.  However,  the  iterative  method  can  still  be  used  as  a  good  check  of  the  accuracy  of  the 
probabilities.  If  the  probability  vector  is  iterated  on  and  DIFF  is  very  small  (say  <  10“ ®),  then  the  iterative 
method  has  converged  and  you  have  very  accurate  results.  When  Po  is  small,  the  recursive  formula  will 
not  yield  very  accurate  results  but  will  still  produce  a  good  starting  vector  for  the  iterative  formula.  When 
Po  <  10-3°,  the  approximation  Po  =  10-3°  is  used  in  order  to  get  some  rough  starting  vector. 

The  system  Mx / D/c  has  an  interesting  property  that  has  been  used  to  speed  up  the  calculations.  If 
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N  is  the  greatest  common  divisor  between  the  different  group  sizes  and  the  number  of  servers,  then  only 
every  Ntk  probability  will  be  non-zero.  Thus  the  program  has  to  only  calculate  every  Ntk  probability. 

Method  for  Finding  the  Roots 

Since  the  basic  root-finding  method  is  discussed  in  Chaudhry  et  al.  (1988),  it  is  briefly  given  here. 

The  roots  of  (1)  that  are  lying  on  the  real  axis  inside  the  unit  circle  are  found  by  using  a  combination 
of  the  bisection  and  Muller’s  method  for  finding  roots.  Then  the  complex  roots  in  the  upper  half  plain  of 
the  unit  circle  are  found  by  using  a  special  algorithm.  In  order  to  use  this  algorithm,  equation  (1)  has  to 
be  modified.  Multiplying  the  R.H.S.  of  (1)  by  e2T*n,  writing  z  =  re**  and  taking  logarithms,  we  can  rewrite 
(l)  as  the  two  equations 

dnr=  h(^(A(z)- i))  (9) 

and 

c9  +  2irn  =  I -  1)^  (10) 

where  i?(u;)  and  I(w)  denote  the  real  and  imaginary  parts  of  w  respectively.  The  algorithm  is  then  for 
j  =  1,2,...,^-*  (where  m  is  the  number  of  real  roots): 

(i)  fix  9  at  9 a  =  Oj-i  +  6  where  90  =  0  and  6  is  an  increment  initially  defined  as  x/(c  +  1)  and  later 
(for  j  >  2)  adjusted  to  minify.!  -  9j_2)/2.0,n/(c  +  l)). 

(ii)  solve  (9)  for  r6a  by  using  a  combination  of  the  bisection  and  Muller’s  methods. 

(iii)  evaluate  the  corresponding  value  of  n^,  the  solution  of  (10)  for  n  using  9a  and  r$a. 

(iv)  fix  9  at  9p  —  6j_x  +  26,  and  obtain  r$9  and  n*9 . 

Defining  f(t)  as  the  fractional  part  of  t,  the  algorithm  now  repeats  with  9a  =  9p  and  8p  =  Op  +  6  if 
/(n#„)  ^  otherwise  there  is  a  root  between  9a  and  9g.  Knowing  there  is  a  root  between  9a 

and  9fi,  we  then  obtain  and  for  9n  =  .5(0a  +  8g),  and  if  >  f(n(o),  set  nta  =  ntj  and 

9a  =  9 otherwise  =  n ^  and  9 g  =  9^.  This  process  of  bisecting  the  sector  (0a,9g)  continues  until 
Qp  ~9a  <  10-4.  Then  the  value  Zy  ■=  r^exp(i0«j)  is  a  root  of  (9)  and  (10).  However,  because  of  the  natural 
logarithmic  function,  this  approximation  needs  to  be  refined  a  bit  more.  Therefore,  this  loot  is  used  as 
starting  value  for  Muller’s  method  using  the  c.e,  which  will  yield  a  more  accurate  root. 
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If  the  angle  being  looked  at  is  greater  than  t  before  all  the  roote  in  the  upper  half  plane  are  found, 
then  the  algorithm  haa  missed  a  root.  In  these  rare  cases,  the  algorithm.is  repeated  using  a  smaller  6  and 
a  smaller  increment  for  the  bisection  method.  This  increases  the  computation  time  but  does  overcome  the 
problem. 


Structure 

SUBROUTINE  MXDC(NNZ,NA,PA,ABAR,C,RHO,MAXIT,DIFM,RTS,P,SIZE) 
Input  Variables: 


NNZ 

Integer 

Number  of  arriving  groups 

NA 

Integer  array(GS) 

Group  Sizes  (indexed  from  i=l,  ...  NNZ  in  order  of  size) 

(eg.  NA(1)=  10,  NA(2)  =  30,  NNZ=2) 

PA 

Real  array(GS) 

Prob(group  size  of  NA(i)  arrives) 

(eg.  PA(1)  =  0.4,  PA(2)  =  0.6) 

Mean  group  size  =  3 
Number  of  Servers  =  c 
Utilization  factor  =  p 
Maximum  number  of  iterations 

Maximum  difference  beween  consecutive  probability  vectors 

Number  of  iterations  performed 
Difference  between  corresponding  probability  vectors 
If  DIFF  =  -1,  then  the  recursive  method  was  used. 

RTS  Complex  axray(RT)  Roots  of  the  c.e. 

P  Real  array(PS)  Prob(number  in  the  system  =  i) 

SIZE  Integer  Size  of  the  probability  vector. 

The  sizes  of  the  arrays  are  specified  in  PARAMETER  statements  in  the  subroutines  and  can  easily  be 
changed.  The  maximum  sizes  of  the  root  vector,  probability  vector  and  number  of  group  sizes  are  in  RS, 


ABAR  Real 

C  Integer 

RHO  Real 

MAXIT  Integer 

DIFM  Real 

Output  Variables: 

MAXIT  Integer 

DIFM  Real 
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PS,  and  GS  respectively.  The  program  has  not  been  thoroughly  tested  for  cases  where  RS  is  larger  than  100. 
GS  can  be  increased  and  should  not  cause  problems.  The  accuracy  of  the  iterative  method  is  dependent  on 
the  sise  of  PS  (see  the  accuracy  section).  The  value  of  3000  for  PS  has  proved  to  be  large  enough  in  most 
cases  tried. 

This  subroutine  calls  upon  the  following  subroutines;  GETRAD,  MULLER,  MULERS,  and  FN.  GE- 
TRAD  is  used  by  the  root-finding  section  of  the  program.  Given  a  particular  angle  and  function,  this  routine 
will  find  the  radius  that  will  give  a  zero  of  the  function.  MULLER  is  a  version  of  Muller’s  method  that  finds 
only  real  roots.  It  is  used  by  GETRAD.  MULERS  is  a  version  of  Muller’s  method  that  finds  complex  roots 
and  is  used  by  the  second  part  of  the  root  finding  section  of  the  program.  FN  computes  values  of  c.e.  and 
real  and  imaginary  parts  of  the  In  of  the  c.e..  This  is  used  by  the  root  finding  section  of  the  program  and  is 
discussed  in  Chaudhry  et  al.  (1988). 
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Accuracy  and  Limitations 

The  accuracy  of  the  root-finding  algorithm  is  measured  by  comparing  the  value  of  the  c.e.  at  the  roots 
to  the  value  of  zc  at  the  roots.  Given  a  root  z,  if  |zc  - exp{6A(A(z)  —  1) } |  <r  |zc|,  then  z  is  a  good  root.  The 
accuracy  of  the  probabilities  is  determined  by  comparing  means  (L)  and  variances  (Var)  as  calculated  by 
using:  the  roots,  the  first  c  probabilities,  and  al!  the  probabilities.  These  formulas  are  discussed  by  Chaudhry 
et  al.  (1988)  and  were  used  to  produce  the  enclosed  table  1.  Another  measure  of  accuracy  is  the  value  of 
the  output  variable  D1FM.  This  output  variable  is  equal  to  the  sum  of  the  absolute  values  of  the  differences 
between  corresponding  probabilities  in  consecutive  probability  vectors.  Thus,  the  smaller  the  value  of  DIFM, 
the  greater  the  convergence  of  the  iterative  method  to  a  steady-state  solution.  In  most  cases  tried,  after  two 
iterations,  the  means  and  the  variances  matched  to  four  or  more  significant  figures,  and  DIFM  <  10~3. 

Two  other  conditions  that  are  required  for  accurate  probabilities  are  that  the  sum  of  the  probabilities 
should  be  very  close  to  one,  and  the  size  of  the  probability  vector  should  not  be  equal  to  the  maximum  vector 
size  PRSIZE.  These  conditions  are  required  for  the  iterative  formula  to  be  valid.  When  p  is  very  high,  the 
size  of  the  probability  vector  can  sometimes  exceed  the  memory  capacity  of  the  computer.  In  these  cases, 
the  prematurely  truncated  vector  will  typically  converge  slowly  to  some  steady-state  solution,  which  may  or 
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may  not  be  close  to  the  actual  solution. 

Time 

The  time  required  to  find  the  roots  ranged  from  only  taking  a  few  seconds  to  taking  a  few  minutes 
(on  an  IBM  PC  AT  with  8087  numeric  co-processor),  depending  on  how  many  roots  had  to  be  found.  The 
time  required  to  calculate  the  probabilities  varied  drastically  depending  on  the  queuing  parameters  and  the 
accuracy  required  For  reasonable  accuracy  (i.e.,  having  the  means  and  the  variances  match  to  within  4 
significant  figures  and  DIF M  <  10-3)  the  time  required  would  be  in  the  order  of  a  few  minutes.  This  time 
would  increase  if  greater  accuracy  was  required,  or  if  very  high  values  for  C  and  p  were  used. 

Precision 

Single  precision  should  be  adequate  for  machines  using  at  least  64  bits  to  represent  real  numbers. 
Otherwise  double  precision  is  required.  In  order  to  change  the  program  to  double  precision,  the  explicit 
REAL  statement  has  to  be  changed  to  DOUBLE  PRECISION,  and  the  functions  REAL,  EXP,  CEXP,  ABS, 
CABS,  AIMAG,  CONJG,  SQRT,  ACOS,  COS,  SIN,  LOG,  replaced  by  DREAL,  DEXP,  CDEXP,  DABS, 
CDABS,  DIMAG,  DCONJ,  DSQRT,  DAGOS,  DCOS,  DSIN,  DLOG.  Also  all  the  constants  throughout  the 
program  have  to  be  changed  to  double  precision  (ie.  1.0  changes  to  1.0D0). 

Example 

A  sample  run  is  given  in  table  1.  All  the  input  variables  and  some  of  the  output  variables  are  shown. 
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Table  1 


Example  of  Input/Output  for  M^/D/10,  with  p-  0.5,  a^-  0.5, 
a2-  0.25,  a3-  0.125,  0.125 


ABAR  -  1.875 
C  -  10 
RHO  -  0.5 
DIFM  -  1.0E-9 


Output  Variables: 

MAXIT  -  5 

DIFM  -  0 . 8702723E-06 
RTS ( 1 )  -  (  1.0000000,  0.0000000) 
RTS ( 2 )  -  (-0.7178427,  0.0000000) 
RTS ( 3 )  -  (  0.5857134,  0.5578024) 
RTS (4)  -  (  0.1643529,  0.7369143) 
RTS (5)  -  (-0.2728908,  0.6757375) 
RTS (6)  -  (-0.5893607,  0.4034265) 
RTS (7)  -  RTS (10)  are  conjugates 
of  RTS (3)  -  RTS (6) 

P(0)  -  0.0630749 

P(l)  -  0.0860505 
P(2)  ■=  0.1021267 


Check  on  Output: 

L  (roots)  -  5.2860534  Var(roots)  -  13.6705954 

L  (first  c  prob . )  =  5.2860540  Var(first  c  prob.)  -  13.6705917 

L  (all  the  prob.)  -  5.2860342  Var(all  the  prob.)  -  13.6704229 

Sum  of  the  Probabilities  =  0.9999989 


Note:  A  user-friendly  software  package  for  the  IBM  PC  family  of  computers 
(or  compatibles)  has  been  developed  by  the  authors.  For  more  information 
about  this  and  other  queuing  packages  please  contact  A&A  publications, 

395  Carrie  Crescent,  Kingston,  Ontario,  K7M-5X7,  Canada. 


P(3)  -  0.1079159 

P(4)  -  0.1161159 

P(5)  -  0.1068249 

P(10)  -  0.0380577 
P(20)  -  0.0008178 
P(30)  -  0 . 8036800E-05 
P(40)  -  0 . 7397169E-07 
P(50)  -  0 . 6576123E-09 
P(60)  -  0 . 5322294E- 11 
P( 70)  -  0 . 2077907E- 13 
P( 72)  -  0 . 1221037E- 13 
SIZE  -  72 


Input  Variables: 


NNZ  -  4 
NA(1)  -  1 
NA(2)  -  2 
NA(3)  =  3 
NA(4)  -  4 


MAXIT  -  5 
PA( 1)  -  0.5 
PA(2)  -  0.25 
PA(3)  -  0.125 
PA(4)  -  0.125 


SUBROUTINE  MXDC (NNZ , NA , PA , ABAR , C , RHO , MAXIT , DI FM , RTS , P , 

6  SIZE) 

C 

C*********************************************************** 

C  This  program  is  designed  to  compute  probabilities  * 

C  of  the  number  in  the  system  for  the  bulk-arrival,  * 

C  multiserver  queue  MAx/D/c.  * 

C***************************************************+******* 

C 

INTEGER  PS.RS.GS 

PARAMETER ( PS-3  000 , RS-5 10 , GS-3 1 ) 

INTEGER  NRR, NCR, I , C , SIZE , NUMIT , MAXIT , IBEG , FLAG 

INTEGER  NNZ, J ,N,M,KSIZE,NA(GS) , IGCD.GCD 

REAL  FI , F2 , FF , X , RHO , ABAR , ANGSP, ANG , PQB(0 : PS) , F , P(0 : PS) 

REAL  PA(GS) ,RAD1 ,RAD2 , ANGO , DIFF , DIFM , CDF , PI (0 : PS) , SUM3 
REAL  RADIUS , ANG1 , ANG2 , ANGLE , TOL, PI 2 , SUM , SUMTOL , RUNSUM 
COMPLEX  RTS (RS) , COEF(RS) ,RT, CF, PC 
PARAMETER (TOL-1 . E- 14) 

C 

C  INITIALIZE  VARIABLES  AND  FIND  THE  REAL  ROOTS 
C 

PI 2-4 . 0*AC0S (0 .0) 

ANG— 2 . 0*ACOS ( 0 . 0 ) / ( C+l ) 

X-0.005 

5  ANGLE-0 . 0 
DIFF— 1.0 
NUMIT-0 
NRR-0 

RADI- -1.0 -TOL 
RAD2— 1 . O-TOL 

10  CALL  GETRAD(GS, RADIUS, ANGLE, RADI ,RAD2,X, 1 , FLAG , F, PA , NA , NNZ , C . RHO 
6  , ABAR . RT , CF) 

IF  ( FLAG . EQ . 1 )  THEN 
NRR-NRR+1 
RTS (NRR) -RAD I US 
RADI— RADIUS+2+X 
GOTO  10 
END  IF 
NRR-NRR+1 

IF(NRR.NE.l)  RTS(NRR)— RTS(l) 

RTS(1)-(1. 0,0.0) 

C 

C  Find  Complex  Roots 


X-X*10 
ANGO-O . 0 
NCR-0 
RADI— 2*X 


mm 


/.-V-  A  V^«r, 


.«VV 


SB9SBH 


o  o  o  non 


t .<■  h“  *»*  1.1  (.<  !>>  i.i  Ii*  ft*  tit 


k * i.t  li • l4  V4» *j‘k  »»  i 


sv-«y«¥ 

v.v 


RAD2-1-T0L 
X-X* . 8 

IF(NRR.EQ.C)  GOTO  90 
ANG1-ANG/2 
start  of  loop 
FI— 0 . 0 

DO  30  ANGLE-ANG1 , PI 2/2 , ANG 

CALL  GETRAD (GS , RADIUS , ANGLE , RADI , RAD 2 , X , 3 , FLAG , FF , PA , NA , 
6  NNZ , C , RHO , ABAR , RT , CF) 

RAD1-RADIUS* . 8 

CALL  FN ( GS , PA , NA , NNZ , C , RHO , ABAR , F2 , ANGLE , RADIUS , 4 , RT , CF ) 
IF(F1.GT.F2)  GOTO  40 
F1-F2 
CONTINUE 
GOTO  60 

ANG1-ANGLE - ANG 
ANG2-ANGLE 

IF(ANG2 -ANG1 . GT . 1 . OE-4)  THEN 
ANGLE- ( ANG1+ANG2 ) /2 . 0 

CALL  GETRAD(GS , RADIUS , ANGLE , RADI , RAD2 , X , 3 , FLAG , FF, PA , NA , 
6  NNZ, C, RHO, ABAR, RT.CF) 

RADI— RADIUS* . 8 

CALL  FN(GS, PA, NA, NNZ, C, RHO, ABAR, FF, ANGLE, RADIUS, 4 ,RT,CF) 
IF(FF.LT.Fl)  THEN 
ANG 2-ANGLE 
F2-FF 
ELSE 

ANG l— ANGLE 
Fl-FF 
END  IF 
GOTO  50 
END  IF 

CALL  MULERS (GS, ANGLE, RADIUS, FF, PA, NA, NNZ, C, RHO, ABAR  RT  CF) 
NCR-NCR+1 
RTS (NCR+NRR)~RT 
IF( (NCR*2+NRR) . LT. C)  THEN 
ANGSP-ANGLE-ANGO 
ANGO-ANGLE 

IF  (ANGSP.LT. ANG)  ANG-ANGSP 
ANGl— ANGLE+ANG/2 
GOTO  20 
END  IF 
end  of  loop 
GOTO  70 
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Missed  some  roots  --  try  again  with  smaller  step  size 


ANG— ANG/ 2 .0 
X-X/50 . 0 
GOTO  5 


FINISHED  ROOT  FINDING  -  COMPUTE  CONJ . 
DO  80  N-l.KCR 

RTS ( NCR+NRR+N ) -CONJ  G ( RTS ( NRR+N ) ) 
CONTINUE 
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c 

C  CALCULATE  THE  PROBABILITIES  P(n) 

C 

C  Incialize  Variables 
C 

90  SUMTOL-1 . 0  -l.E-12 
SIZE-0 
IBEG-0 

DO  100  J-O.PS 
PQB(J)-0.0 
PI(J)-0.0 
P(J)-0.0 
100  CONTINUE 

DO  110  I-1,C 
COEF(I)-0 .0 
110  CONTINUE 
C 

C  Find  the  Greatest  Common  Divisor 
C 

GCD-0 

DO  130  N-NA(l) ,1,-1 
DO  120  J-l.NNZ 

IF( (NA(J)/N)*N.NE.NA(J) )  GOTO  130 
120  CONTINUE 

IF(GCD.EQ.O)  GCD-N 
IF( (C/N)*N . EQ . C)  GOTO  140 
130  CONTINUE 
140  IGCD-N 
C 

C  Calculate  Pi's 
C 

KSIZE-PS 

PI ( 0 ) -EXP ( - C*RHO/ABAR ) 

SUM-PI (0) 

C 

DO  170  N-GCD, (PS -1) ,GCD 
PI (N)-0 . 0 
DO  150  I-l.NNZ 

IF( (N-NA(I) ) . LT. 0)  GOTO  160 
PI(N)-PI(N)+NA(I)*PA(I)*PI(N-NA(I) ) 

150  CONTINUE 

160  PI (N)— PI (N)*C*RHO/(ABAR*N) 

SUM— SUM+PI (N) 

IF(SUM.GE. SUMTOL)  THEN 
SUM-SUM- PI (N) 

KSIZE-N 
GOTO  180 
END  IF 
170  CONTINUE 

KSIZE— N-GCD 
180  DO  190  N-KSIZE , PS 
PI(N)— 0.0 
190  CONTINUE 
C 

C  CALCULATE  P(0) 

C 


ilMS  »■* 


PC-1.0 

DO  230  I  -  2.C 

PC  -  PC*RTS ( I)/(RTS (I ) -1.0) 

230  CONTINUE 

P(0)  -REAL(PC)*C*(1 .O-RHO) 

IF(P(0) .LE.1.0E-30)  THEN 
P(0)-1 . 0E-30 
END  IF 
C 

C  CALCULATE  THE  COEFICIENTS  FOR  THE  FIRST  C  PROBABILITIES 
C 

C0EF(1)  -  - 1 . EO/RTS (1) 

DO  250  I  -  2,  C 
DO  240  J  -  2,  I 

C0EF(I - J+2)  -  C0EF(I-J+2)  -  COEF(I-J+l)/RTS(I) 
240  CONTINUE 


COEF(l)  -  C0EF(1) 
CONTINUE 


-  1. EO/RTS (I) 


C  CALCULATE  P(l)  --->  P(C-l) 

C 

SUH-P(O) 

IF  (C.EQ.l)  GOTO  270 
DO  260  I  -  IGCD ,  C-l.IGCD 
P(I)-REAL(COEF(I) )*P(0) 

IF(P(I) .LT.0.0. OR . P ( I ) .GT.1.0)  THEN 
SIZE— I 
P(l)-0.0 
GOTO  300 

END  IF 
SUM— SUM+P ( I ) 

I F ( SUM . GE . SUMTOL) THEN 
SIZE-I 
GOTO  300 
END  IF 
260  CONTINUE 
C 

C  CALCULATE  P(C) , P(C+1) ,P(C+IGCD)  -  - ->  P(SIZE) 
C 

270  P(C)-=P(0)*EXP(C*RHO/ABAR) -SUM 
IF(P(C) .LT.0.0)  P(C)-0.0 
SUM-SUM+P(C) 

P(C+1)— (P(l) -PI ( 1)*SUM)/PI (0) 
IF(IGCD.EQ.l)  THEN 
IBEG-2 
ELSE 

IBEG-ICCD 
END  IF 

CDF— SUM+P (C+l) 

DO  290  N-IBEG, (PS-C) , IGCD 
IF(N- IGCD . LT . KSIZE)  THEN 
I-C+ICCD 
ELSE 

I-N+C- KSIZE 
ENDIF 
SUM3-0 . 0 


I 
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DO  280  M-I , C+N- 1 , IGCD 

SUM3-SUM3+P(M)*PI(N-M+C) 

280  CONTINUE 

P(C+N)— (P(N) -PI (N)*SUM-SUM3)/PI (0) 
IF(P(C+N) . LT . 0 . 0) P(C+N)-0 . 0 
CDF— CDF+P(C+N) 

IF( CDF. GE.SUMTOL) THEN 
SIZE— C+N 
GOTO  300 
ENDIF 

290  CONTINUE 
SIZE-PS 

PREPARE  FOR  ITERATION  LOOP 

300  SUM-0 . 0 

DO  310  N-0 , SIZE 
PQB(N)— P(N) 

SUM— SUM+PQB(N) 

310  CONTINUE 

IF(SUM.LT. 0.9999)  THEN 
DO  320  N-0 .SIZE 

PQB(N)— PQB(N) /SUM 
P(N)— PQB(N) 

320  CONTINUE 

ENDIF 


START  OF  ITERATION  LOOP 

DO  390  NUMIT-1,  MAXIT 
SUM-0 . 0 

DO  330  I-O.C-l 

SUM-SUM+PQB ( I ) 

330  CONTINUE 

P(0)— PI (0)*(SUM+PQB (C) ) 

RUNSUM-P(O) 

DO  360  N-IGCD , PS , IGCD 
SUM3-0 . 0 

DO  340  J-=MINO(KSIZE,N)  ,0,-1 

IF((N-J+C) .GT. (SIZE-1))  GOTO  350 
SUM3-SUM3+PQB(N-J+C)*PI(J) 

340  CONTINUE 

350  P(N)=PI(N)*SUM+SUM3 

IF( (RUNSUM+P(N) ) .GT. SUMTOL)  GOTO  370 
RUNSUM-RUNSUM+P ( N ) 

IF(P(N) .LT.l.E-14. AND. RUNSUM.GT. 0.9999)  GOTO  370 
360  CONTINUE 

370  SIZE— N- IGCD+1 

DIFF-0.0 

DO  380  N-0, SIZE- 1, IGCD 

DIFF-ABS ( P(N) -PQB(N) )+DIFF 
PQB(N)— P(N) 

380  CONTINUE 

IF(DIFF. LE.DIFM)  GOTO  400 
390  CONTINUE 


ooo  uooo  U  o  u 


NUMIT-NUMIT - 1 


RETURN  TO  MAIN  PROGRAM 

*00  DIFM-DIFF 

MAXIT-NUMIT 

RETURN 

END 


SUBROUTINE  GETRAD (GS , RADIUS , ANGLE , RADI , RAD 2 , X , N , FLAG , FF , PA , NA , 
6  NNZ , C , RHO , ABAR , RT , CF) 


INTEGER  GS 

INTEGER  FLAG, N , NA(GS ), NNZ, C 

REAL  ANGLE , RADIUS , RADI , RAD2 ,X, FF, FI , RHO, ABAR , PA(GS ) 

COMPLEX  RT , CF 
C 

FLAG-1 
RADIUS-RAD 1 

CALL  FN(GS, PA, NA, NNZ, C, RHO, ABAR, FI, ANGLE, RADIUS, N.RT.CF) 

10  RADIUS-RADIUS+X 

IF(RADIUS .GT.RAD2)  GOTO  20 

CALL  FN(GS, PA, NA, NNZ, C, RHO, ABAR, FF, ANGLE, RADIUS, N.RT.CF) 
IF(FF*F1 . LT . 0 . 0)  THEN 
RADIUS-RADIUS -X/2 . 0 
GOTO  30 
END  IF 
Fl-FF 
GOTO  10 
20  FLAG-0 

RADIUS-1. 0-1. 0E- 10 

3C  CALL  MULLER(GS, ANGLE, RADIUS, FI, FF.N.X, PA, NA, NNZ, C, RHO, ABAR, RT.CF) 

40  RETURN 

END 
C 
C 

C  . 

c 

SUBROUTINE  MULERS (GS , ANGLE , RADIUS , FF , PA , NA , NNZ , C , RHO , ABAR , RT , 

6  CFRTS) 

C 

C  . 

c 

C  THIS  ROUTINE  IS  THE  COMPLEX -VALUED  FORM  OF  MULLERS 
C  METHOD  USED  TO  FIND  ROOTS  OF  THE  CHAR.  EQN . 

C 

INTEGER  GS 

INTEGER  NA(GS) , NNZ , C 

COMPLEX  Cl , H , Z , DELFPR , CFRTS , FRTPRV , FRTDEF , LAMBDA 

COMPLEX  DELF , DFPRLM , NUM , G , SQR , DEN , RT 

REAL  RADIUS , PI2 , ANGLE , TOL1 , TOL3 , FF , RHO , ABAR , PA ( GS ) 


PARAMETER (Cl- (0 . 0,1-0) .TOL3-1 . 0E- 9 , T0L1-1 . OE- 16) 

PI2— 4*ACOS(0.0) 

KCOUNT— 1 
KOUNT-O 

KCOUNT— KCOUNT  + 1 
IF(KCOUNT . EQ . 7)  GOTO  40 
MAXIT-25 

Z-RADIUS*CEXP ( C I * ANGLE ) 

H- . 00001E0*10 . EO**KCOUNT 
RT-Z+H 

CALL  FN (GS, PA, NA.NNZ , C ,RHO , ABAR , FF , ANGLE, RADIUS, 2 ,RT,CFRTS) 

DELFPR-CFRTS 

RT-Z-H 

CALL  FN(GS, PA, NA.NNZ, C.RHO, ABAR, FF, ANGLE, RADIUS, 2, RT.CFRTS) 

FRTPRV-CFRTS 

RT-Z 

CALL  FN(GS , PA , NA , NNZ , C , RHO , ABAR , FF , ANGLE , RADIUS , 2 , RT , CFRTS ) 

FRTDEF-CFRTS 

DELFPR-FRTPRV - DELFPR 

LAMBDA— 0.5  EO 

DELF— FRTDEF- FRTPRV 

DFPRLM-DELFPR* LAMBDA 

NUM-- FRTDEF* (1 . 0+LAMBDA)*2 . 0 

G— (1 . 0+LAMBDA*2 . 0) *DELF-LAMBDA*DFPRLM 

SQR— G*G+2 . 0*NUM* LAMBDA* ( DELF - DFPRLM) 

IF( (REAL(SQR) . LT . 0 . 0) .AND. (ANGLE. GT. (PI2/2 . EO-O . 001E0) ) ) 

6  SQR-(0. 0,0.0) 

SQR-CSQRT ( SQR) 

DEN-G+SQR 

IF(ANGLE . GT . (PI2/2 . 0-0 . 001E0) )THEN 
I F ( ( REAL ( G ) *REAL ( SQR ) ) . LT . 0 . 0 ) D EN-G - SQR 
GOTO  20 
END  IF 

I F ( REAL ( G ) *REAL ( SQR ) +AIMAG ( G ) *AIMAG ( SQR ) . LT . 0 . 0 ) D  EN-G - SQR 

IF(CABS(DEN) . LT . T0L3 ) DEN-1 . EO 

LAMBDA-NUM/DEN 

FRTPRV-FRTDEF 

DELFPR-DELF 

H-H*LAMBDA 

RT-RT+H 

IF(KOUNT.GT.MAXIT)  GOTO  40 
KOUNT— KOUNT + 1 
I F ( CABS ( RT ) .CT.1.2E0)  THEN 

ANGLE-ACOS ( REAL (RT) /CAB S(RT) ) 

RT— . 99999999*CEXP(CI*ANGLE) 

ENDIF 

CALL  FN(GS , PA , NA , NNZ , C , RHO , ABAR , FF , ANGLE , RADIUS , 2 , RT , CFRTS ) 
FRTDEF-CFRTS 

IF (CABS (H) . LT . T0L1*CABS (RT) )  GOTO  40 

IF( CABS ( CFRTS ) .LT. T0L1)  GOTO  40 

IF(CABS ( FRTDEF) . LT . 10. 0*CABS( FRTPRV) )  GOTO  10 

H-H/2.0 

LAMBDA-LAMBDA/2 . 0 
RT-RT - H 
GOTO  30 


40  RADIUS-CABS (RT) 

ANGLE-ACOS ( REAL(RT ) /RADIUS ) 

RETURN 

END 

C 

C 

C  . 

c 

SUBROUTINE  MULLER (GS , ANGLE , RADIUS , FI , FF , N , X , PA , NA , NNZ , C , RHO , ABAR , 
6  RT.CF) 

C 

C  . 

c 

C  THIS  ROUTINE  IS  THE  REAL-VALUED  FORM  OF 

C  MULLERS  METHOD 

C 

INTEGER  GS 

REAL  FRTPRV , DELFPR , FRTDEF , DELF , DFPRLM , G , SQR , H , DEN 
REAL  LAMBDA ,NUM, FI, RHO, ABAR, PA(GS) , RADIUS .ANGLE , FF , X 
INTEGER  KOUNT,MAXIT,N,NA(GS) ,NNZ,C 
COMPLEX  RT , CF 
C 

KOUNT-O 
MAXIT-25 
H— X*0 . 5E0 
FRTPRV-F1 
DELFPR-F1 - FF 

CALL  FN(GS , PA, NA, NNZ, C, RHO, ABAR, FF, ANGLE, RADIUS, N, RT.CF) 
FRTDEF-FF 
LAMBDA- -0.5E0 
10  DELF-FRTDEF- FRTPRV 

DFPRLM-DELFPR* LAMBDA 

NUM—  FRTDEF* ( 1 . 0+LAMBDA) *2 . 0 

G- ( 1 . 0+LAMBDA*2 . 0 ) *DELF - LAMBDA*DFPRLM 

SQR-G*G+2 . 0*NUM*LAMBDA*( DELF -DFPRLM) 

I F ( SQR . LT . 0 . 0 ) THEN 
LAMBDA-G/ ( 2 . 0*LAMBDA* ( DFPRLM - DELF) ) 

ELSE 

SQR-SQRT(SQR) 

DEN-G+SQR 

IF(G*SQR . LT . 0 . 0) DEN-G - SQR 
IF(ABS(DEN) .LT.l.E-7)  DEN-1.0 
LAMBDA-NUM/DEN 
ENDIF 

FRTPRV— FRTDEF 
DELFPR-DELF 
H— H*LAMBDA 
RADIUS— RADIUS+H 
IF(KOUNT.GT.MAXIT)  GOTO  30 
C 

20  KOUNT— KOUNT+1 

CALL  FN(GS , PA , NA , NNZ , C , RHO , ABAR , FF , ANGLE , RADIUS , N , RT , CF) 

FRTDEF-FF 

C  CHECK  FOR  CONVERGENCE . 

IF(ABS(H) .LT. l.E-14)  GOTO  30 
IF(ABS( FRTDEF) .LT.l.E-16)  GOTO  30 
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C  CHECK  FOR  DIVERGENCE . 

IF(ABS (FRTDEF) . LT . 10 . 0*ABS (FRTPRV) )  GOTO  10 
H-H/2 . 0 

LAMBDA-LAMBDA/ 2 . G 
RADIUS-RADIUS -H 
GOTO  20 
C 

30  RETURN 
END 


SUBROUTINE  FN(GS , PA,NA, NNZ , C , RHO , ABAR , F , ANGLE , RADIUS , N , RT , CF) 


INTEGER  GS 

INTEGER  NA(GS) ,NNZ , I , C ,N 

REAL  SUM, PA(GS ) ,F, RADIUS , ABAR, ANGLE, RHO, RN,TOL,PI2 
COMPLEX  CSUM.RT.CF 

TOL-l.E-14 
PI2— 4*ACOS(0.0) 

IF(ABS (RADIUS) .LT.1.0E-4)  RADIUS-1 . OE-4 

SUM-0.0 

CSUM-<0. 0,0.0) 

GOTO  (100,200,300,400)  N 

real  version  of  char,  equation 

DO  110  I-l.NNZ 

SUM-SUM+PAC I ) * (RADIUS**NA (I) ) 

CONTINUE 

F— RADIUS**C - EXP( (C*RHO/ABAR)*(SUM- 1 . 0) ) 

CF-F 

GOTO  500 

complex  version  of  char,  equation 

IF(CABS (RT) . GT .1.01)  THEN 
F— . 1E0 
GOTO  500 
END  IF 

DO  210  I-l.NNZ 

CSUM— CSUM+PA( I ) *(RT**NA( I ) ) 

CONTINUE 

CF-RT**C - C  EXP ( ( C*RHO/ABAR ) * ( C  SUM -1.0)) 

F— CABS (CF) 

GOTO  500 

modified  char,  equation 
DO  310  I— 1 , NNZ 

SUM— SUM+PA(I)*(RADIUS**NA(I) ) *COS (NA( I ) *ANCLE) 
CONTINUE 
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F- ( C*RHO/ABAR) * ( SUM- 1 . 0 ) - C*LOG (RADIUS ) 

GOTO  500 

m  value 

DO  410  I-l.NNZ 

SUM-SUM+PA (I)*( RADIUS  **NA < I ) ) *S IN ( NA ( I ) * ANGLE ) 
CONTINUE 

RN-((C*RH0/ABAR)*SUM-C*ANGLE)/PI2 

I-RN+TOL 

RN-RN-I 

F-ABS(RN) 

RETURN 

END 


RK? 

Bra 


non  n  a  o  o  o  o 


PROGRAM  INPUT 


INTEGER  PS.RS.GS 

PARAMETER( PS-3000, RS-510 , GS-31) 

INTEGER  C,C2,I,NNZ,NA(GS) .MAXIT, SIZE, U 

REAL  RHO , RRHO , Al , A2 , A3 , SUM3 , SUM4 , LI , L2 , L3 , VI , V2 , V3 

REAL  PA(GS) , P(0: PS) , ABAR , SUM1 , SUM2 

REAL  RADIUS, ANGLE ,RN,F,DIFM, SUM 

COMPLEX  RT(RS) ,R, CF 

WRITE(*,*)  '  WELCOME  TO  THE  QUEUEING  SYSTEM  MXDC' 
WRITE(*,*) 

WRITE(* , *) 

WRITE(* , *)  'PLEASE  ENTER  THE  #  OF  DIFFERENT  GROUP  SIZES' 
READ(*,*)  NNZ 
WRITE(* , *) 

WRITE(* , *) '  ENTER  THE  GROUP  SIZE  DISTRIBUTION:' 

ABAR-0 . 0 
DO  100  I-l.NNZ 

READ ( * , * )  NA ( I ) , PA ( I ) 

ABAR-NA(I)*PA(I)  +  ABAR 
100  CONTINUE 

WRITE(* , *) 

WRITE(* , *)  '  ABAR  -  ' , ABAR 
WRITE(*,*) 

WRITE(* , *)  '  ENTER  RHO  ' 

READ ( * , * )  RHO 
WRITE(* , *) 

WRITE(* , *)  '  ENTER  #  OF  SERVERS  (C) ' 

READ(* , *)  C 
WRITE(* , *) 

URITE(* , *)  '  ENTER  DIFM' 

READ ( * , * )  DIFM 
WRITE(*,*) 

WRITE(* , *)  '  ENTER  THE  MAX  #  OF  ITERATIONS  (MAXIT)’ 

READ(* , *)  MAXIT 

WRITE(*,*) 

WRITE(* , *)  '  CALCULATING' 

U-6 


CALL  MXDC (NNZ, NA, PA, ABAR, C, RHO, MAXIT, DIFM ,RT,P, SIZE) 
C 
C 
C 
C 

WRITE(U , *)  '  THE  ROOTS  ARE:’ 

URITE(U , *) 

WRITE(U , *) 

DO  800  I-l.C 


RADIUS-CABS  (R) 

ANGLE-ACOS (REAL(R) /RADIUS ) 

CALL  FN(GS, PA, NA.NNZ.C.RHO, ABAR.RN, ANGLE .RADIUS, 4 ,R,CF) 
CALL  FN(GS, PA, NA,NNZ,C,RHO,ABAR,F, ANGLE,  RADIUS, 2, R.CF) 
WRITE(U,*)  I,'  ROOT  -  ' ,RT(I) 

VRITE(U,*)  '  RN  -  ' ,RN 
WRITE(U,*)  '  FN  -  ' ,F 

WRITE(U,*) 

800  CONTINUE 


W 

sSsft 


RRHO-l-RHO 
C2— C*C 
Al— 0.0 
A2-0.0 
A3-0.0 

DO  1000  I-l.NNZ 

A1-NA(I)*PA(I)+A1 
IF(NA(I) .GT.l)  THEN 

A2-NA(I)*(NA(I)-1)*PA(I)  +  A2 
IF(NA(I) .GT.2)  THEN 

A3-NA(I)*(NA<I)-1)*(NA(I)-2)*PA(I)+A3 

ENDIF 

ENDIF 

CONTINUE 


SUM1-0 . 0 
SUM2-0 . 0 
SUM3-P(0)*C2 
SUM4— P(0)*C2*C 

DO  1100  I-l.C 

IF(I.NE.l)  THEN 

SUMl-REAL(l/(l-RT(I)))  +  SUM1 
SUM2-REAL(RT(I)/((1-RT(I))*(1-RT(I))))+SUM2 
ENDIF 

IF(I.NE.C)  THEN 

SUM3-(C2-I*I)*P(I)  +  SUM 3 
ENDIF 

SUM4-(C2*C- 1*1*1 ) *P(  I ) +SUM4 
CONTINUE 


LI— SUM1+1/ ( 2*RRHO)*(RHO*A2/Al+l - C*RRHO*RRHO) 

VI— 4*RH0*A3*RRH0/A1  +  3*RHO*A2/Al*RHO*A2/Al 
Vl-Vl+6*RHO*A2/Al*(C*RHO*RHO- 2*C*RHO+C-RHO+2 )  +1+2*RH0 
Vl-(Vi+6*C*RHO*RRHO*RRHO-C2*RRHO*RRHO*RRHO*RRHO) 

Vl-Vl/( 12*RRH0*RRH0)  -  SUM2 

L2— (C*RHO*( 1+A2/A1)  -  C2*RRH0*RRH0  +SUM3 ) / ( 2*C*RRH0) 

V2-- 3*L2*(C2*RRHO*RRHO+C*( 1 - 2*RH0) - C*RH0*A2/Al ) 
V2-V2+3*C*RHO*A2/Al*(-C*RHO+C+l)  +  C*RH0*A3/Al 
V2-V2+C2*C*RHO*RHO*RHO+3*C*C*RHO*RRHO*(C+l)  -C*(C*C-RH0)+SUM4 
V2— V2/( 3*C*RRH0)  +L2-L2*L2 
L3-0 . 0 
V3-0.0 


w 


WWW) 


.  .W  AW.Ii'A'A'A*. 


vi  v.twv-'V'  rJ>.  v*  .v'rj'Tj'  ~j>  urw*j> 


e 


R 

R 


V 

p 

s 


DO  1200  I-0.SIZE-1 
L3— I*P(I)+L3 
V3-I*I*P(I)+V3 
1200  CONTINUE 

V3-V3-L3*L3 
C 
C 

WRITE(U,*)  ' 

WRITE(U,*) 

WRITE(U,*)  ' 

WRITE(U,*)  ' 

WRITE(U, *)  ' 

WRITE(U,*) 

WRITE(U,  *)  ' 

WRITE(U , *)  ' 

WRITE(U , *)  ' 

WRITE(U ,  *) 

WRITE (U , *)  ' 

WRITE(U ,  *)  ' 

WRITE(U , *)  ' 

WRITE(U , *) 

WRITE(U,*)  'MAXIT  -  ' .MAXIT 
WRITE(U, *)  'DIFF  -  ' , DIFM 
WRITE(U , *)  'SIZE-  SIZE-1 
WRITE(U, *) 

WRITE(U,*)  '  THE  PROBABILITIES  ARE:' 
WRITE<U,*) 

SUM-0.0 

DO  1400  1-0 , SIZE- 1 
WRITE(U,*)  '  P(' ,1,')  -  ' ,P(I) 

SUM— SUM+P(I) 

1400  CONTINUE 
WRITE(U , *) 

WRITE(U , *)  '  THE  SUM  OF  THE  PROBS  -  '.SUM 

1401  CLOSE(ll) 

END 


CHECK  ON  OUTPUT' 

USING  ROOTS:' 

MEAN-  ' , LI 
VARIANCE-  ' ,V1 

USING  THE  FIRST  C  PROBS:' 

MEAN-  ’ , L2 
VARIANCE-  ' ,V2 

USING  ALL  THE  PROBS:' 

MEAN-  ' , L3 
VARIANCE-  ' ,V3 


If 

:♦ 


distribution  list 


Office  of  Naval  Research 
800  North  Quincy  Street 
Arlington,  VA  22217 

Attention:  Scientific  Officer, 

Statistics  and  Probability 
Mathematical  Sciences  Division 

ONR  Resident  Representative 
Joseph  Henry  Building,  Room  623 
2100  Pennsylvania  Avenue,  N.W. 

Washington,  DC  20037 

Director,  Naval  Research  Laboratory 
Washington,  DC  20375 

Attention:  Code  2627 

Defense  Technical  Information  Center 
Building  5,  Cameron  Staeion 
Alexandria,  VA  22314 

C.  M.  Harris 

GMU  Office  of  Research 


5$ 


m 


'•  P.WPiI'i'Vi 


WC'V 


